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Abstract 

We have developed a physics-based land-surface temperature (LST) algorithm for simultaneously retrieving 
surface band-averaged emissivities and temperatures from day/night pairs of MODIS (Moderate Resolution 
Imaging Spectroradiometer) data in seven thermal infrared bands. The set of 14 nonlinear equations in the 
algorithm is solved with the statistical regression method and the least-squares fit method. This new LST 
algorithm was tested with simulated MODIS data for 80 sets of band-averaged emissivities calculated from 
published spectral data of terrestrial materials in wide ranges of atmospheric and surface temperature conditions. 
Comprehensive sensitivity and error analysis has been made to evaluate the performance of the new LST 
algorithm and its dependence on variations in surface emissivity and temperature, upon atmospheric conditions, 
as well as the noise-equivalent temperature difference (NEAT) and calibration accuracy specifications of the 
MODIS instrument. In cases with a systematic calibration error of 0.5%, the standard deviations of errors in 
retrieved surface daytime and nighttime temperatures fall between 0.4-0.5°K over a wide range of surface 
temperatures for mid-latitude summer conditions. The standard deviations of errors in retrieved emissivities in 
bands 31 and 32 (in the 10-1 2.5 Jim IR spectral window region) are 0.009, and the maximum error in retrieved 
LST values falls between 2-3 °K. Several issues related to the day/night LST algorithm (uncertainties in the 
day/night registration and in surface emissivity changes caused by dew occurrence, and the cloud cover) have 
been investigated. The LST algorithms have been validated with MODIS Airborne Simulator (MAS) dada and 
ground-based measurement data in two field campaigns conducted in Railroad Valley play a, NV in 1995 and 

1996. The MODIS LST version 1 software has been delivered. 



1. Introduction 


Land-surface temperature (LST) is one of the key parameters in the physics of land-surface processes on regional 
and global scales, combining the results of all surface-atmosphere interactions and energy fluxes between the 
atmosphere and the ground [Mannstein, 1987; Sellers et al., 1988], Therefore it is required for a wide variety of 
climatic, hydrological, ecological and biogeochemical studies [Schmugge and Andre, 1991; Running et al., 
1994], In order to understand the Earth system better on the global scale, the Earth Observing System (EOS) will 
provide surface kinetic temperatures at specified accuracies of 0.3 °K for oceans and 1 °K for land. The 
international Tropical Ocean Global Atmosphere (TOGA) program has specified that sea surface temperature 
(SST) should be accurate to 0.3 °K for global numerical models of climate. 

During the past decade, significant progress has been made in estimation of land-surface emissivity and 
temperature from airborne thermal infrared data. Kahle et al. [1980] developed a technique to estimate the 
surface temperature based on an assumed constant emissivity in one channel and previously determined 
atmospheric parameters. This temperature was then used to estimate the emissivity in other channels [Kahle, 
1986]. Other techniques, such as thermal log residuals and alpha residuals, have been recently developed to 
extract emissivity information from multispectral thermal infrared data [Hook et al., 1992], 

A variety of split-window methods have been developed to retrieve sea-surface temperature and land-surface 
temperature from NOAA AVHRR data. The split-window LST method corrects the atmospheric effects based on 
the differential absorption in adjacent infrared bands [Price, 1984; Becker, 1987; Wan and Dozier, 1989; Becker 
and Li, 1990; Sobrino et al., 1991; Vital, 1991; Kerr et al., 1992; Ottle and Stoll, 1993; Prata, 1994; Wan and 
Dozier, 1996], A major problem in using split-window LST methods is that we need to know the surface 
emissivities in the bands better than 0.01. It seems possible to have such knowledge of the emissivities for 
certain types of land covers, such as lake surfaces, snow/ice, dense evergreen canopies, and some soils. For land 
covers with variable emissivities, especially in semi arid and arid areas, it is almost impossible to estimate two 
band-averaged emissivities to such accuracy. Therefore, it is necessary to develop new algorithms to retrieve 
LST without prior knowledge of surface emissivities. In addition, surface emissivity is also needed both to 
calculate up-welling thermal infrared radiation, and for environmental monitoring and geological mapping 
[Kahle, 1980], 
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Li and Becker [1993] proposed a method to estimate both land-surface emissivity and LST using pairs of 
day/night co-registered AVHRR images. They used a temperature-independent spectral index (TISI) in thermal 
infrared bands and assumed knowledge of surface TIR BRDF (thermal infrared Bidirectional Reflectance 
Distribution Function) and atmospheric profiles. Such combined a priori knowledge and information are not 
readily available in most situations. Wan and Li [1996] developed a MODIS day/night LST method that uses 
day/night pairs of TIR data in seven MODIS bands for simultaneously retrieving surface temperatures and band- 
averaged emissivities without knowing atmospheric temperature and water vapor profiles to high accuracy. 

MODIS is an EOS instrument that will serve as the keystone [Salomonson et al., 1989] for global studies of 
atmosphere [King et al., 1992], land [Running et al., 1994], and ocean processes. It scans ±55° from nadir in 36 
bands, with bands 1-19 and band 26 in the visible and near infrared range, and the remaining bands in the thermal 
infrared from 3 to 15|im. The specifications of MODIS bands are shown in Table I. Figure 1 shows the total 
transmission and transmission functions corresponding to water vapor (H 2 O) band absorption and continuum 
absorption, uniformly mixed gases (labeled by CO 2 ) and ozone absorptions. These transmission calculations 
were done for a typical scenario, using an updated atmospheric radiative transfer model, MODTRAN3 [Berk et 
al., 1989], the specific conditions are mid-latitude summer, clear-sky (23-km visibility), with total column water 
vapor of 2.9cm and line-of-sight viewing of 45°. The transmission functions corresponding to molecular 
scattering, aerosol scattering and absorption are also shown in this figure. By examining the atmospheric 
transmission functions, we can have some general ideas for applications of the MODIS bands. Those bands in 
transparent atmospheric windows are designed for remote sensing of surface properties. Other bands are mainly 
for atmospheric studies. The exact location and bandwidth of MODIS bands are selected to meet the 
requirements from atmospheric, ocean and land sciences. MODIS will provide images of daylight reflection and 
day/night emission of the Earth, repeating global coverage every 1-2 days. It uses 12 bits for quantization in all 
bands. The thermal infrared bands have an IFOV (instantaneous field-of-view) of approximately 1 km at nadir. 
MODIS will view cold space and a full-aperture blackbody before and after viewing the Earth scene in order to 
achieve calibration accuracy of better than 1 % absolute for thermal infrared bands. MODIS is particularly useful 
because of its global coverage, radiometric resolution and dynamic ranges, and accurate calibration in multiple 
thermal infrared bands designed for retrievals of SST, LST and atmospheric properties. Specifically, all 
atmospheric channels of MODIS will be used to retrieve atmospheric temperature and water vapor profiles 
[Smith et al., 1985], band 26 will detect cirrus clouds, and thermal infrared bands 20, 22, 23, 29, 31-33 will 
correct for atmospheric effects and retrieve surface emissivity and temperature. Multiple bands in the mid- 
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infrared range will provide, for the first time, corrections for solar radiation in daytime LST estimations using 
mid-infrared data. Because of its multiple bands in the mid-infrared range and in the 8- 14 pm window, MODIS 
provides an unprecedented opportunity to develop a physics-based algorithm to simultaneously retrieve surface 
emissivity and temperature. In Section 2, we present the theoretical basis of the new LST algorithm. Section 3 
describes numerical methods used in the algorithm. Section 4 gives simulation results of the new LST algorithm 
and results from sensitivity and error analysis. In Section 5, we investigate several issues related to the day/night 
LST algorithm (uncertainties in day/night registration and in surface emissivity changes caused by dew 
occurrence, and the cloud cover). In Section 6, we present validation results of the LST algorithms with MODIS 
Airborne Simulator (MAS) dada and ground-based measurement data in two field campaigns conducted in 
Railroad Valley playa, NV in 1995 and 1996. Finally, the status of MODIS LST software is described in Section 
7. 

2. Theoretical Basis of the New LST Algorithm 

For land covers with variable and unknown emissivities, there is greatly insufficient mathematically under- 
determined information to retrieve surface temperature and band-averaged emissivities from a one-time 
measurement of N thermal infrared channels. This is so, even when atmospheric temperature and humidity 
profiles are known exactly, there are still N + 1 unknowns (N band emissivities plus surface temperature). 
Therefore, we will consider using multi-temporal and multi-channel data. 

2.1 A Physics-Based Day/Night LST Model 

In clear-sky conditions, the spectral infrared radiance L(K, p) at the top of the atmosphere is composed of surface 
thermal emittance, thermal path radiance L a (k, p), path radiance resulting from scattering of solar radiation 
L S (X, p, Po, <t>o), solar beam and downward solar diffuse radiation and atmospheric thermal radiation reflected by 
the surface, 

L(k, p) = 1 1 (A,, p) e(X, p) B(A, T s ) + L a (X, p) + L s (k, p, po, <j>o) + 1 2 (X, p, pp) po Eo(A.)/ r (p; po, <t>o) 

2n 1 

+ J JV/r(|LP', p)L</(A., -p', 0') + r 4 (A, p)L,(A.,-p', <)>')] dp'd<t>' , (1) 

0 0 

where p is cosine of the viewing zenith angle, e(X, p) is the surface spectral emissivity, B(X,T 5 ) is the radiance 
emitted by a blackbody at surface temperature T s , Eq(K) is the spectral solar irradiance incident on the top of the 
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atmosphere (normal to the beam), po is cosine of the solar zenith angle, is the relative azimuth between the 
viewing direction and the solar beam direction, / r (p; p', §') is the BRDF function, Lj(k y -p\ §') is the downward 
solar diffuse radiance, L t (X, -p', <|f) is the atmospheric downward thermal radiance, their incident direction is 
represented by -\i' and <))', and ?/(),/ = 1, . . , 4 are transmission functions for the corresponding terms. 

The wavelength, X, in eq. (1) is the wavelength center of a narrow wavelength interval because there is no way to 
measure the exact monochromatic signal as a continuous function of wavelength by satellite sensors. Equation 
(1) is the generalized form used in the 8-14 pm thermal infrared range [Wan and Dozier, 1990] into a wider 
wavelength range of 3- 14pm. It requires complete calculations of the atmospheric radiative transfer to determine 
the values of all terms on the right-hand side. After the zenith and azimuth dependent radiance at user-defined 
levels from the Earth’s surface to the top of the atmosphere (TO A) is provided by accurate atmospheric radiative 
transfer simulations, the TOA radiance can be represented by its components in form (1). Its special form has 
been used for a long time in many atmospheric radiation models including LOWTRAN [Kneizys et al., 1983], 
MODTRAN [Berk et al., 1987], and MOSART models [Comette et al., 1994]. In the special form, 
ti(X, p) = 1 1 ( X , p) and t^(X y p) = t \ ( X , p) are assumed. 

In order to retrieve surface emissivity and temperature from eq. (1), we need to use suitable TIR bands. 
According to the MODIS band specifications in Table I and the atmospheric transmission in Fig. 1, bands 20, 22, 
and 23 are in the transparent atmospheric window in the 3. 5-4. 2 pm medium wavelength range, bands 29-32 are 
in the 8-13 pm atmospheric window range, while band 33 is just on the edge of this atmospheric window. Band 
30 is strongly affected by ozone absorption, so using this band does not help to retrieve surface temperature. As 
shown in Fig. 1, major absorbers in bands 20, 22, and 23 are CO 2 , N 2 , and water vapor. Major absorbers in bands 
29, and 31-33 are water vapor and CO 2 . The transmission corresponding to aerosol scattering and absorption in 
these bands is in the range of 0.95-0.98. So using average aerosol distribution in atmospheric radiative transfer is 
usually good enough unless volcano eruptions strongly change the aerosol distribution. Since CO 2 and O 2 
mixing ratios are almost constant, their densities are determined by atmospheric pressure and temperature. Water 
vapor is the most variable absorber in the Earth’s atmosphere. Therefore, if we know atmospheric water vapor 
and temperature profiles, we can calculate all atmospheric terms in the above equation to a quite high accuracy, 
which is limited mainly by the accuracy of the coefficients of the water vapor continuum and band absorptions. 
The MODIS sounding channels can be used to retrieve atmospheric temperature and water vapor profiles [Smith 
et al., 1985; Menzel and Purdom, 1994]. But retrieving atmospheric profiles needs the knowledge of surface 



6 


emissivity in order to separate the surface contribution from the sounding data. Therefore, the retrieved profiles 
might not be accurate in areas where surface emissivities are highly variable such as in semi-arid and arid areas. 
Although the absolute values of the retrieved profiles are not accurate, the shapes of the atmospheric temperature 
and water vapor profiles may be reasonably well obtained. Radiative transfer simulations show that the radiance 
at the top of the atmosphere, in MODIS TIR bands 20, 22, 23, 29, 31-33, is almost not affected by changing 
atmospheric temperature and water vapor profiles at levels above 9km elevation because the values of water 
vapor content at these levels are very low. If the shapes of temperature and water vapor profiles in the lower 
troposphere can be well retrieved from the MODIS sounding data, we can use two variables to describe the 
atmospheric variations. One is the amount of shift in the temperature profile up to 9km elevation. Another is the 
scale factor for the water vapor profile so that we can determine the column water vapor with the shape and the 
scale factor. Then we use the atmospheric temperature at the surface level, T a , as the representative variable of 
the tropospheric temperature profile. Similarly, we can use the column water vapor (cwv) as the representative 
for the water vapor profile. Alternatively, we can consider it as the first order of approximation to describe the 
atmospheric condition by using these two variables. 

In order to make practical use of multi-temporal and multi-channel data, we need to simplify eq. (1) by using 
some realistic assumptions about the surface optical properties. We assume: 1) The surface emissivity changes 
with vegetation coverage and surface moisture content. However, it does not significantly change in several days 
unless rain and/or snow occurs during the short period of time - particularly for bare soils in arid and semi-arid 
environments (for which the surface of the ground is normally dry) [Kerr et al., 1992]. 2) There are quite strong 
spectral variations in surface reflectance for most terrestrial materials in the 3. 5-4. 2pm wavelength range 
[Salisbury and D’Aria, 1994], but their BRDF anisotropic factor in this wavelength range has very small 
variations in the order of 2% [Snyder and Wan, 1996(b); Snyder and Wan, 1996(c)]. So it seems appropriate to 
assume that a single BRDF anisotropic factor can be used for the surface-reflected solar beam in MODIS bands 
20, 22 and 23 located in this wavelength range. This anisotropic factor is defined by the ratio of the surface- 
reflected solar beam at the view direction of the MODIS sensor to the radiance that would have resulted if the 
surface reflected isotropically (such a surface is called Lambertian surface), 

7i/ r (p;po,^o) 

a = , (2) 

r 

where r is reflectance of the assumed Lambertian surface. 3) Atmospheric radiative transfer simulations show 
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that in clear-sky conditions the surface-reflected diffuse solar irradiance term is much smaller than the surface- 
reflected solar beam term in the thermal infrared range, and the surface-reflected atmospheric downward thermal 
irradiance term is smaller than surface thermal emission. So the Lambertian approximation of the surface 
reflection does not introduce significant error in the 3-14|im thermal infrared region. Then we can replace the 
BRDF function f r ( p, | <JT) in eq. (1) with r/n and link it to the surface emissivity e by r = 1 - e according to 
Kirchhoff’s law. 

It is important to point out that in eq. (1) we separate the surface-reflected solar beam term from its irradiance 
term (the integral of the downward solar diffuse radiance) because changing solar zenith angle has different 
effects on these two terms. As solar zenith angle increases, the solar beam at the surface level decreases, but the 
downward solar diffuse irradiance may increase in some situations. If the solar beam is included in the total solar 
irradiance incident on a surface and surface reflectance (also called as hemispherical reflectance or albedo in the 
visible and near-infrared range) is defined as the ratio of the total solar radiance reflected from the surface to the 
total solar irradiance, the surface reflectance will be dominated by the BRDF of the solar beam and therefore the 
reflectance significantly depends on solar zenith angle [Schaaf and Strahler, 1993]. After the solar beam is 
separated from the total downward solar irradiance, we can use the BRDF anisotropic factor to calculate the 
surface-reflected solar beam and use the surface reflectance to calculate the surface-reflected solar downward 
irradiance. In this way, the solar angle and viewing angle dependences in the surface reflectance will be smaller 
so that we can assume the surface as a Lambertian surface. 

Based on above assumptions, we have developed the following physics-based day/night LST model from eq. (1). 
The radiance measured in MODIS band j can be expressed as 

L(j) = ti(J)e(j)Bj(T s ) + L a (J) + L s (j ) + 1 ^ Ll1 [r 2 0) « Ho £oO) + t 3 (J)E d (j) + U(j) E,(j )] , (3) 

K 

where all terms are band-averaged, e(j) is the band emissivity which will be given in eq. (5), similarly for B j(T s ) y 
E a (j), L s {j), and Eq(j ). Ej(j) and E t (j) are the band-averaged solar diffuse irradiance and atmospheric 
downward thermal irradiance at surface. And tj(j) , I = 1,.., 4 are the band effective transmission functions 
weighted by the band response function, the corresponding radiance, and irradiance terms. Note that we have 
neglected the in-band spectral variation of the surface emissivity in reducing eq. (1) into eq. (3), and have omitted 
symbols of view angle and solar angle for most terms in the above equation. On the right-hand side of this 
equation, £(/'), a, and B j(T s ) depend on surface properties and conditions, all other terms depend on atmospheric 
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water vapor and temperature profiles, solar angle and viewing angle. These terms can be given by numerical 
simulations of atmospheric radiative transfer. The spectral response functions measured from the Engineering 
Model of the MODIS instrument have been used as weights in calculations of band averages of these terms. 

If we use two measurements (day and night) in N MODIS TIR bands, we have 2N observations. The number of 
unknown variables are N band emissivities, daytime surface temperature T s -day* nighttime surface temperature 
T s -night* f° ur atmospheric variables (T^ and cwv at two times), and the anisotropic factor a, totalling N + 7. The 
number of observations must be equal to or larger than the number of unknowns, 

2N >N + 7 . (4) 

So N > 7 . Note that it is necessary to apply independent shapes of atmospheric temperature and water vapor 
profiles for daytime and nighttime so that temporal variations and temperature inversion (more often at night) 
could be considered in the LST retrieval. For the MODIS LST algorithm, these seven bands are MODIS bands 
20, 22, 23, 29, 31-33. According to the experience from the Engineering Model of the MODIS instrument, the 
NEAT in band 33 may be reduced from 0.25 °K to 0.12°K, and it appears possible to achieve the goal for 
absolute calibration accuracy, 0.5-0.75%, for these seven TIR bands. It seems that we can get unique solutions 
for the above 14 unknowns using 14 observations. But it is actually not true because: 1) the atmospheric profile 
is a continuous function of height and there are only a finite number of MODIS sounding bands so that the 
atmospheric temperature and water vapor profiles can be retrieved only at a finite number of levels, 2) there are 
always uncertainties in the retrieved atmospheric profiles and even in their shapes, 3) there is always instrument 
noise in the measurement data, 4) there are uncertainties in the atmospheric optical properties including water 
vapor absorption coefficients which we used in the development of LST algorithms. Therefore all we can do is to 
use a best combination of available bands and use an appropriate method to determine the best estimates of the 
unknown variables. We also need to use enough a priori knowledge and constraints of the atmosphere and the 
surface as "virtual measurements" to make the retrieval problem well posed [Rodgers, 1976]. The advantage of 
using daytime data in MODIS bands 20 and 22-23 is that solar radiation can be used as TIR source in the medium 
wavelength range so that the day/night LST model is essentially an active method to get the information of 
surface reflectance. Combining with the nighttime data in these three bands and day/night data in the other four 
MODIS bands makes it possible to simultaneously retrieve surface emissivity and temperature. The advantage of 
including four atmospheric variables ( T a and cwv in daytime and nighttime) is that atmospheric variations are 
considered in the retrieval procedure so that uncertainties in the initial atmospheric conditions could be reduced 
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and a better self-consistent solution of the surface emissivity and temperature could be reached. 
2.2 Band-Averaged Emissivities of Land- Surface Materials 


The band-averaged emissivity is defined as 


£(/) = 


Xj.u 

J *F(X) e(X) dX 
hhk 


j 'V(X)dX 

h-L 


(5) 


where X F(^) is the spectral response function of band j, Xj i and Xjjj are its lower and upper boundaries. By 
using MODIS spectral response functions, band-averaged emissivities can be calculated from published spectral 
reflectance data of 80 pure terrestrial materials [Salisbury and D’Aria, 1992 and 1994]. This spectral data base 
includes igneous, metamorphic, and sedimentary rocks, varnished rock surfaces, lichen-covered sandstone, soil 
samples, green foliage, senescent foliage, water, ice, and water surfaces with suspended quartz sediment or oil 
slicks. The sample names and numbers are listed in Table II. The calculated band emissivities in MODIS bands 
20, 22, 23, 29, 31-33 are shown in Fig. 2. The sample number corresponds to the sample name and the type of 
material in Table II. As shown in this figure, there are very strong variations in the band emissivities for rock and 
sand samples, and for some soil samples and senescent vegetation foliages. For example, the emissivity of sands 
in MODIS band 20 could be as low as 0.55. However, the band emissivities in MODIS bands 31-33 are larger 
than 0.8 for all samples in the spectral reflectance database. For water, ice, and green vegetation leaves, there are 
small emissivity contrasts among these seven bands and their band emissivities vary in small ranges. 


2.3 Atmospheric Radiative Transfer Simulations 


The accuracy of atmospheric radiative transfer numerical models depends on numerical methods to solve the 
radiative transfer equation and both our knowledge of the atmosphere and its inherent optical properties [Goody 
and Yung, 1989]. The many methods available to solve the atmospheric radiative transfer problem have proved 
their fundamental theory and mathematically interesting because there are important applications in neutron 
diffusion theory, astrophysics, and earth sciences. For example, there are a variety of methods based on two- 
stream approximations [Meador and Weaver, 1980; Tanre et al., 1990], 4-stream approximations [Cuzzi et al., 
1982], and others such as delta-M method [Wiscombe, 1977], adding/doubling method [Wiscombe, 1976], 
discrete ordinate method [Stamnes and Conklin, 1984], and Monte Carlo simulation method [Adams and 
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Kattawar, 1978]. We developed a radiative transfer model, which provides accurate matrix solutions of the 
azimuth-dependent scalar radiative transfer equation for a vertical inhomogeneous, multi-layer atmosphere by 
using the adding/doubling method for the development of algorithms to estimate column ozone and LST in 
clear-sky conditions [Smith et al., 1992; Wan and Dozier, 1989; Wan and Dozier, 1996]. Results from this model 
match those from Stamnes and Conklin’s discrete-ordinates [Stamnes and Conklin, 1984] calculations to four 
decimal places. Atmospheric radiative models based on the adding/doubling method have advantages in easy 
implementation of surface interfaces, such as the air-water interface and interfaces for specular reflectance or 
BRDF reflectance, and in efficiently getting solutions for multiple boundary conditions. 

The continuous update of the LOWTRAN code [Kneizys et al., 1983; Kneizys et al., 1988] and MODTRAN code 
[Berk et al., 1987; Berk et al., 1989] developed by the U.S. Air Force Geophysics Laboratory over the past two 
decades represents a significant progress in improving our knowledge of optical properties of the earth’s 
atmosphere. A two-stream approximation with multiple scattering parameterization [Kneizys et al., 1988; Isaacs 
et al., 1987] is used in the LOWTRAN model. The new versions of MODTRAN code have the option to use the 
discrete ordinate method. 

In the TIR range, LOWTRAN7 [Kneizys et al., 1988] and MODTRAN [Berk et al., 1989] give transmission 
functions of each molecule at a wavenumber interval of 5cm -1 and 1cm -1 , respectively. The LOWTRAN7 band 
model was based on a least squares fitting to FASCODE [Clough et al., 1986] calculations of various 
transmittances, but MODTRAN is based on a separate band model, derived directly from the HITRAN [Rothman 
et al., 1992] data base, by-passing FASCODE. The MODTRAN band model has been validated against 
spectrally degraded FASCODE calculations for both radiance and transmittance. The agreement is usually 
within a few percent RMS and seldom exceeds 5 % [Wang et al., 1996]. Note that non-monochromatic 
transmission values cause a violation of the Lambert-Bouguer-Beer law because of the complexity of molecular 
band absorption. This is true even for a narrow wavenumber interval of 1 cm -1 . One solution to this problem is 
to expand radiative transmission functions calculated from LOWTRAN or MODTRAN by using “exponential- 
sum fitting” [Wiscombe and Evans, 1977]. Following this technique the monochromatic radiative transfer model 
[Chandrasekhar, 1960] is applied separately to each term in the exponential-sum expansion, and the results 
summed. After convoluting these results with the spectral response function, we can get band-averaged 
atmospheric terms in eq. (3). A radiative transfer code using the “exponential-sum fitting” formulation has been 
implemented successfully on workstations and massive parallel computers, and the simulation results have been 
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used in the development of a generalized split-window LST algorithm [Wan and Dozier, 1996]. Typically, the 
total number of the cross-product terms in the exponential-sum fitting formulation is over 1,000 so it is very 
computationally time consuming. The advantage in using the exponential-sum fitting formulation is that we 
obtain more accurate results. For example, the resulted three effective transmission functions for the viewing 
path in eq. ( 1 ) may be different by a few percent (£3 >t\ and cfj) because of selective, wavelength- 
dependent molecular band absorption. This occurs even though these transmission functions are defined for a 
same optical path from the target to the top of the atmosphere in the viewing direction of the sensor. The 
correlated-k distribution method [Lacis and Oinas, 1991; Fu and Liou, 1992] is an alternative to the exponential- 
sum-fitting method. Usually the number of terms used in the correlated k-distribution method is smaller so that it 
is computationally efficient and it gives accurate results. The multiple scattering algorithm in the MODTRAN 
code is being upgraded to include a “correlated-k” absorption characterization. 

As shown in Fig. 1, the atmospheric transmission in the 8- 13 Jim window, where three MODIS bands are located 
for the purpose of remote measurements of surface temperature, strongly depends on water vapor absorption 
including band absorption and continuum absorption. A review for measurements of water vapor absorption in 
the 8-13 |Lim atmospheric window reveals a considerable variation in its magnitude over the past 20 years [Grant, 
1990]. The accuracy of water vapor continuum absorption in five of the measurements reviewed is 
approximately 10%, adequate experimental measurements are lacking at temperatures below 280 °K. There is no 
accepted theory for the continuum absorption. Recent theoretical studies [Ma and Tipping, 1992; Ma and 
Tipping, 1994] on water vapor continuum absorption have led to significant progress in understanding the 
physical mechanisms and the temperature dependence of the continuum absorption. But it is still premature to 
theoretically determine the magnitude and the temperature dependence of the water vapor continuum absorption 
coefficients. Thus, modelers must rely on empirical formulations [Robert et al., 1976; Clough et ah, 1989] based 
on laboratory measurements [Burch and Alt, 1984], Atmospheric conditions, especially cold temperatures and/or 
high humidities, are difficult if not impossible to reproduce in the laboratory. This is particularly true in the vital 
area of continuum absorption. Studies at relative humidities over 70% are a persistent problem. This is the 
threshold for condensation on hygroscopic dust particles and therefore for fogging of optical elements. 
Furthermore, laboratory spectroscopists have almost reached an impasse in the area of line wings and the 
continuum that prevents progress in line-by-line modeling [Varanasi, 1988]. In the past several years, the water 
vapor absorption has been compared and validated with High-Resolution Interferometer Sounder (HIS) spectral 
radiance data involving vertical path measurements from an aircraft, as well as from the ground [Clough et al.. 
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1989; Smith et al., 1993], and long-path atmospheric transmission measurements [Theriault et al., 1994]. The 
empirical continuum formulation used in the FASCODE code [Clough et al., 1986], as well as in LOWTRAN and 
MODTRAN codes has been changed several times in the past decade. In the LOWTRAN7 code [Kneizys et al., 
1988] and earlier versions of the MODTRAN code [Berk et al., 1987], the magnitude of the water vapor self 
continuum absorption coefficient in the 8.5-1 3 (xm window is smaller by approximately 20 % compared to that 
used in the LOWTRAN 6 code [Kneizys et al., 1983]. Clough [1995] made a new correction to the the water 
vapor continuum based on the measurement of the downwelling radiance at Kavieng, New Guinea by Westwater 
et al. [1994], and the measurements by Revercomb and colleagues at the University of Wisconsin. This new 
continuum formulation has been implemented in version 3 of the MODTRAN code in 1994. Because the 
uncertainty in the water vapor continuum absorption coefficients may be larger than a few percent, it is not 
critical to neglect the small differences among transmission functions t\,t 3 and *4 in eq. ( 1 ) before the accuracy 
of the water vapor continuum absorption coefficients is significantly improved. So we accept the approximations 
*3 = 0 and *4 = t\ in eq. ( 1 ) that are made in the MODTRAN3 code up to its version 1.3 in the time being. The 
effect of these approximations should be considered along with uncertainties due to other sources when actual 
data are used to retrieve surface emissivity and temperature. Note that even we assumed that in each narrow 
spectral interval of 1 cm -1 or 5 cm -1 the transmission functions for the viewing path are equal in eq. ( 1 ), the band 
effective transmission functions t 3 (j) and t 4 (j) may still be different from t\(j) in eq. (3). Keeping in mind all 
the problems raised above for radiative transfer models, in this study we use version 1 .3 of the MODTRAN3 code 
[Berk et al., 1989] to calculate all atmospheric and solar terms in eq. (3). The discrete ordinate option with eight 
streams is used in MODTRAN3 calculations so that the effect of multiple scattering due to background aerosols 
is considered in the calculations of the path radiances and the downwelling irradiances. 

The path radiance resulting from scattering of the solar radiation in eq. (3), L s (j), does depend on the relative 
azimuth between viewing direction and the solar beam direction. This dependence is determined by the aerosol 
loading, its size distribution, type and scattering phase function. The aerosol information and properties are not 
readily available in most situations. As shown in Fig. 1 , the total aerosol effect on the transmission function in 
the thermal infrared range is small in normal clear-sky conditions. Radiative transfer simulations indicate that 
the value of L s is only several hundredths of the direct solar beam value at the surface level, and that the azimuth 
dependence in L s is less than 10%. So it is appropriate to neglect this azimuth dependence and to use the 
azimuth-averaged value of L s in the new LST algorithm. 
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2.4 Variations of Atmospheric Conditions 

It is important that a practical LST algorithm should accommodate atmospheric variations in a range that is wide 
enough to cover all possible real situations. For LST retrieval, we only consider atmospheric variations in clear- 
sky conditions. In the thermal infrared range, the most important atmospheric variables are atmospheric water 
vapor and atmospheric temperature profiles. Atmospheric absorption and thermal emission occur mainly in the 
lower troposphere. We assume that the MODIS product of the atmospheric temperature and water vapor profiles 
retrieved from MODIS sounding channel represents the shapes of the profiles well but maybe does not represent 
their absolute values because of the difficulties in decoupling the atmosphere-land interaction. Atmospheric 
temperature and water vapor at any level will be interpolated from their values retrieved at fixed levels. Given 
the shapes of temperature and water vapor profiles, we can use only two variables to describe variations of the 
clear-sky atmospheric condition: a shift of the temperature profile below elevation 9km, and a scaling factor for 
the water vapor profile. The column water vapor can be determined by the shape and the scaling factor. In order 
to build a data base for the atmospheric and solar terms in eq. (3), we will select 24 basic atmospheric profiles 
considering different shapes of temperature and water vapor profiles, and the range of air surface temperatures in 
different regions and seasons. Some basic atmospheric profiles include temperature inversion layers. Then we 
add more variations to each of these basic atmospheric profiles in the following ways: 1) A 57 is added to the 
atmospheric temperature profile at all levels between surface and elevation 9km, 57 varies from -10°K to +20 °K 
in steps of 2°K. The modified atmospheric temperature at the surface level, T a , will be used as representative of 
the entire atmospheric temperature profile. 2) The atmospheric water vapor profile at all levels between surface 
and elevation 9km is scaled in steps of 10% so that the column water vapor varies from 10% to 120% of the basic 
value. 

2.5 Variations of the Land-Surface Temperature 

In the simulation study of the new LST algorithm, we consider LST variations in a wide range of LST variations. 
The daytime surface temperature varies from atmospheric surface temperature T a ^ ay to T a -^ ay + 24 °K in steps 
of 6°K, and the nighttime surface temperature varies from T a ^ n i g } lt - 13.5°K to T a - n i g h t +4.5°K in steps of 


4.5 °K. 
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3 . Numerical Methods Used In The New LST Algorithm 

3.1 Look-Up Table Method 

In order to save computational time on numerical simulations of atmospheric radiative transfer for calculating the 
atmospheric and solar terms in eq. (3), look-up tables will be used in the new LST algorithm. In this way, we 
only need to make a complete series of radiative transfer simulations once to build these look-up tables. Because 
multi-dimension interpolations are involved in our look-up table method, linear interpolation is most efficient. 
This requires smaller intervals (or steps) for these look-up tables. For example, the step for the atmospheric 
temperature variation is 2°K, the step for atmospheric column water vapor is 10% of the average value, the step 
for solar zenith angle and viewing zenith angle is 10° for angles smaller than 30° and 5 ° for larger angles. The 
upper limits for solar and viewing zenith angles are 75° and 65°, respectively. Similarly, a look-up table in a step 
of 0.1 °K is also built for the band-averaged Planck functions in the temperature range 200-400 °K, which is the 
maximum dynamic range for the seven MODIS bands used in this LST algorithm. It is required that errors 
caused by look-up tables and interpolation methods should be smaller than NEAT. If this resolution scheme is 
used to build a look-up table for the three solar terms in eq. (3) for 24 basic atmospheric profiles, the total size of 
the look-up table is approximately 60MB. If use of small look-up tables is a higher priority, we can use 3-point 
interpolation method so that six zenith angles are enough for viewing and solar angles. However, 3-point 
interpolation takes much more computational time than linear interpolation. The sizes of look-up tables for the 
other three atmospheric terms (transmission, thermal infrared path radiance, and downward thermal infrared 
irradiance) are much smaller. 

3.2 Approaches to Solve the Retrieval Problem 

We have developed two approaches to solve the LST retrieval problem. The first one uses statistical regression 
method, and the second one uses other numerical methods to solve the set of nonlinear eq. (3). 

In a linear approximation of eq. (3) in the proximities of reference values of surface temperature and band 
emissivities, the left-hand side reduces to the band brightness temperature and the right-hand side reduces to 
surface temperature and band emissivities. Combining 14 equations together, the solution for surface 
temperature and band emissivities should be a linear combination of the band brightness temperatures, each of 
which corresponds to one of the 14 observations. Its mathematical form is 



15 


14 

x i - X W i.j yj + W <\0 , ( 6 ) 

7 = 1 

where x is a vector of the 14 variables including surface temperatures and band emissivities, yj is the band 
brightness temperature for observation j, and i = 1, ..., 14 and j = 1, ..., 14 are coefficients. And w; o is the 
coefficient for the offset term. We can determine these coefficients in two steps. In step 1 , we construct a large 
set of simulated observation values in wide ranges of atmospheric and surface conditions. In step 2, we make a 
statistical regression analysis using the band brightness temperatures associated with these simulated band 
radiance values as independent variables and using the given surface band emissivities and temperatures, and 
atmospheric parameters as dependent variables. The output of this regression analysis will be the coefficients in 
eq. (6). The process of statistical regression analysis takes much computational time. But it needs only to be 
done once. The values of jc f provided by this approach are the best estimates of these unknown variables in the 
statistical sense. 

If we have better information on the shapes of the atmospheric temperature and water vapor profiles for the time 
which makes it possible to have a clear-sky day/night pairs of MODIS data, we can use other methods to 
numerically solve the set of nonlinear eq. (3). We tried the Quasi-Newton method [Dennis and Schnabel, 1993] 
and the Least-Squares Fit (% 2 fit) method [Bevington 1969]. As Rodgers [1976] pointed out, retrieval problems 
in remote sensing are generally nonlinear. The main sources of the nonlinearity in eq. (3) are: 1) temperature 
dependence of the atmospheric transmission, 2) the dependence of transmission on absorber concentration, 3) 
temperature dependence of the Planck function, 4) wavelength dependence of the Planck function across a 
spectral band, 5) wavelength dependence of the Planck function between spectral bands, and 6) nonlinear 
constraints. 

The initial values of the 14 unknown variables are given in their constrained ranges based on reasonable guesses 
or statistical analysis. The Quasi-Newton method is slightly more computationally efficient. These two methods 
give similar results in cases not including noise. It is well known [Bevington, 1969; Dennis and Schnabel, 1983] 
that global convergence to right solutions is not guaranteed for nonlinear problems, especially when noise is 
included. The % 2 fit method is selected in the new LST algorithm because it is more stable in our simulation 
studies. We are only interested in real situations where there is noise in remote measured data caused by the 
intrinsic instrument sources and actual turbulence and fine structure in the atmosphere. 
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A measure of the goodness of y} fit is defined by [Bevington, 1 969] 


X 2 = Z {^-[Lj-L{j)Y) , 
7 = 1 


14 


( 7 ) 


where Lj is the scaled band radiance observation value, j = 1, 7 for daytime, 7 = 8 , 14 for nighttime. L(/) is the 
scaled band radiance function in eq. (3), which depends on unknowns jq , i = 1 , 14. We use the values of band- 
averaged Planck functions at a reference temperature, 300 °K, to scale the band radiance in corresponding bands 
so that the scaled differential radiance may be comparable. The term Oj is the uncertainty in observation value 
Lj. In cases without noise, Cj is identically equal to 1. However, for cases which include noise NEAT, Oj will be 


Cj = Lj 


n j t^neq (./ ) 

T b U) 


( 8 ) 


based on the following approximation for the band-averaged Planck function 

Lj ~Cj Tfr nj (j ) (9) 

where A T neq {j) is the NEAT value in band j, and T^ij) is the brightness temperature corresponding to band 
radiance Lj. In the temperature range 240-400 °K, regression analysis gives the best fitting values for tij, they are 
12.91, 12.25, 11.98, 6.00, 4.70, 4.11, and 3.74 for MODIS bands 20, 22, 23, 29, 31-33. Note that this 
approximation is used only in calculation of Oj , which determines the weight in eq. (7). The effect of errors 
caused by this approximation on solutions is negligible. 

One of the difficulties in the fit processing is that there may be more than one local minimum for within a 
reasonable range of values for variable jq, particularly in cases including noise. Therefore the final solution may 
depend on their initial values. We use two different ways to make the initialization. In the first way, we use a 
dozen sets of initial values that are spread over preassigned ranges all from minimums to maximums to get 
different solutions and select the solution associated with the minimum % value. In noisy situations, this 
selected solution may not be the best one we searched for. An alternative way is to use the estimates provided by 
the statistical regression method as the initial values. We use the second way in our LST algorithm. Typically, 
the x 2 fit method takes 3-4 iterations to reach the final solution. 
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4. Sensitivity and Error Analysis 

Using look-up tables, we can quickly construct 14 band radiance values (seven values for daytime and seven 
values for nighttime) at the top of the atmosphere for any given conditions of surface band emissivities and 
BRDF anisotropic factor, daytime and nighttime surface temperatures, daytime and nighttime atmospheric 
surface temperatures and column water vapor values, solar angle and viewing angle. Then we can use these 14 
radiance values as simulated MODIS observations to retrieve the given surface and atmospheric variables. We 
can then construct a sensitivity and error analysis, presented in the following sub-sections. 

4.1 Errors Due to Look-up Table and Interpolation Methods 

In the first numerical simulation experiment, we do not include any noise in the data construction in order to test 
the numerical method to solve the nonlinear problem and to evaluate the errors caused by the use of look-up 
tables and interpolation methods. We use the temperature and water vapor shapes in the “average” mid-latitude 
summer atmospheric profiles (model 2 in MODTRAN) and set the daytime and nighttime atmospheric surface 
temperatures at 298.2 °K and 290.2 °K. The column water vapor is set at 2.6cm for both daytime and nighttime 
for simplicity. In real applications, we use independent variables for the column water vapor in daytime and 
nighttime. We set the anisotropic factor as 1 , solar zenith angle at 45 °, viewing angle at nadir for daytime and 
nighttime, five different daytime surface temperatures ranging from 298.2 °K to 322.2 °K, and five different 
nighttime surface temperatures ranging from 276.7 °K to 294.7 °K. There are 25 cases of different daytime and 
nighttime surface temperatures for each sample of 80 surface materials. The band emissivities of these 80 
terrestrial material samples cover the range from 0.55 to almost unity. The standard deviations of errors in 
retrieved surface temperatures are 0.27 and 0.21 °K for daytime and nighttime, the standard deviations of errors in 
retrieved emissivities are in 0.005-0.008 for bands 1 to 6, and 0.012 for the last band because of the low 
transmission of MODIS band 33 in the atmospheric condition. The standard deviations of errors in retrieved 
BRDF anisotropic factor, atmospheric temperatures, and column water vapor are 0.08, 0. 10-0.1 5 °K, and 0.06cm, 
respectively. These numbers indicate that look-up tables are appropriate and the fit method works well. 

4.2 Sensitivities to the Uncertainties in Atmospheric Profile Shapes 

In the second simulation experiment, we set the NEAT values for the seven bands at 0.05, 0.07, 0.07, 0.05, 0.05, 
0.05, and 0.12 K, set 0.5% as the systematic calibration error for all bands, and keep all other parameters as in 
the first experiment. In our simulation, NEAT is treated as a random noise. We consider four different 
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atmospheric conditions in mid-latitude summer, one is the “average” condition used in MODTRAN code (model 
2), two (labeled by A 109 and A117) are selected from the satellite TOVS Initial Guess Retrieval (TIGR) 
atmospheric profile database [Moine et al„ 1987], and the fourth (labeled by “average-4K”) is the variant of the 
“average” one by shifting -4°K on the temperature profile but keeping its water vapor profile unchanged. As 
shown in Fig. 3, three of them have almost the same air temperature at the surface level, but they have very 
different shapes in the temperature and water vapor profiles. The temperature discrepancy between the 
“average” profile and profile A109 may be as large as 10°K at elevations near 2km and between 6-10km. The 
difference in water vapor profiles in atmospheric conditions of “average”, A109, and A1 15 may be 20% to 50% 
or even larger. We established separate data bases of the atmospheric terms in eq. (3) through atmospheric 
radiative transfer simulations for these different atmospheric conditions. The separate data bases will be used to 
calculate the daytime and nighttime band radiances in the seven MODIS bands in wide ranges of surface 
temperature for 80 surface samples. These calculated band radiances are then used as simulated observations. 
The coefficients in eq. (6) were obtained by statistical regression analysis of the observations simulated for the 
“average” atmospheric condition. We suppose that there is enough information available for the average 
atmospheric condition, but there is no information available on the shapes of the atmospheric profile for 
atmospheric conditions likes those in A 109 and A 115. In the statistical approach, we apply the same set of 
regression coefficients to the four sets of simulated observations data for retrieving surface temperatures and 
emissivities. In the y} fit approach, these surface temperatures and band emissivities retrieved by the regression 
approach are used as initial values for further iterative processing. The standard deviations of errors in surface 
temperatures and band emissivities retrieved by the statistical regression method are given in the first part of 
Table III, and those retrieved by using the y 2 fit method are given in the second part. Comparing the results from 
the statistical approach and y} fit approach for the “average” atmospheric condition indicated that the y 2 fit 
method gives significant improvements on retrieved surface temperatures and band emissivities. This is because 
we know the shapes of the atmospheric temperature and water vapor profiles well enough to select the right set of 
the regression coefficient and the right parts from the look-up tables for the atmospheric and solar terms in eq. (3). 
But for atmospheric conditions A109 and A1 15, we do not have the information for making the right selections. 
So the results retrieved from the y 2 fit approach are worse than those from the statistical approach. However, for 
the case with the shifted “average” temperature profile, the standard deviations of errors in surface temperatures 
retrieved by using the y 2 approach is reduced by a factor of 2, and the accuracies of retrieved band emissivities 
are improved by almost 50% because the shapes of atmospheric temperature and water vapor profiles in this case 
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are as same as those in the “average” profile. From this experiment we gained the following insights: the 
statistical method is less accurate but is also less sensitive to uncertainties in the atmospheric profile shapes; and 
the fit method may be more accurate but is more sensitive to uncertainties in the profile shapes. In the 
following part of this paper, we assume that the information of the profile shapes is available so that it is 
appropriate to pursue the % fit approach. 

43 Sensitivities to the Three Assumptions of Surface Optical Properties 

In the first test, Cl , of the third simulation experiment, we set the NEAT values for the 7 bands at 0.05, 0.07, 0.07, 
0.05, 0.05, 0.05, and 0.12°K, set 0.5% as the systematic calibration error for all bands, and kept all other 
parameters as in the first experiment. The errors in surface temperatures retrieved by the % 2 fit method for a total 
of 2,000 different cases are shown in Fig. 4A. The errors in retrieved band emissivities in MODIS bands 31 and 
32 are shown in Fig. 4B. The standard deviations of errors in retrieved surface daytime and nighttime 
temperatures are in range of 0.4-0.5°K, and the standard deviations of errors in band-averaged emissivities in 
MODIS bands 31 and 32 are 0.009 over a wide range of surface temperatures in the mid-latitude summer 
atmospheric condition. We can see the effect of the 0.5% systematic calibration error in Fig. 4A. This forces the 
retrieved temperature to shift to the positive direction by approximately 0.2 °K. The histograms of errors in 
retrieved surface temperatures and emissivities in bands 31 and 32 for a total of 2,000 cases are shown in Fig. 5 A 
and 5B. We also applied the day/night LST method to “mixed samples” that contain a mixture of two samples of 
the 80 terrestrial materials at different surface temperatures. Similar results have been obtained as long as band 
emissivities of the mixed sample are calculated from the components with their proportions as weights and its 
effective surface temperature is calculated from the total thermal radiation in MODIS band 31. The effective 
surface temperature of a mixed sample does depend on band number, but this dependence is very small (about the 
instrument noise level) in normal surface conditions. We do not consider forest fires in the LST processing 
because the MODIS TIR bands may be easily saturated by subpixel fires at a small size and there is no sufficient 
knowledge on the optical properties of fires and smoke for making atmospheric corrections. In test Cl of this 
experiment, the maximum error in the retrieved LST is 3.2 °K for only one case (for the Indian grass sample at 
temperature 322.2 °K). Without considering this extreme case, the maximum error in retrieved LST will be 
2.2 °K. Note that we simulated the surface temperature variation in a very wide range. For each sample, the 
surface temperature varies in a range of 24 °K in daytime, and in a range of 18°K in nighttime. These ranges are 
too wide for some land covers in real situations. For example, the temperature of snow cover and ice could not be 
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above a few degrees C. This is so, even considering the possibility of some small proportion of other land covers 
mixed in the scene, because the temperature of water surface and healthy vegetation is not likely to be warmer 
than the air surface temperature by 24 °K. Higher solar elevation and surface temperature in a normal range are 
favorable conditions for land covers with low reflectances. In such situation, the solar beam can be effectively 
used as an active TIR source for remote measurements of the surface reflectance in the medium wavelength range 
so that the band emissivities in MODIS bands 31 and 32 can be also retrieved at an accuracy better than 0.01. 
Then these retrieved emissivities can be used in the generalized split-window LST algorithm [Wan and Dozier, 

1996] to quickly retrieve LST of same area for a period of one or more weeks depending on season and weather 
conditions. 

Now we check whether it is possible to relax the three assumptions of surface optical properties made in section 
2.1. The first row in Table IV gives standard deviations of the surface temperature and emissivities retrieved in 
test Cl of this experiment. In tests C2 and C3, we introduce some variations for the nighttime surface band 
emissivities to simulate its possible change with surface moisture content. In test C2, the nighttime band 
emissivities increase by 0.01 and they are only limited by its maximum value 1. In test C3, the emissivity 
increment depends on its value, a lower band emissivity could increase more. This may be the case for sands, its 
emissivity in MODIS band 20 is 0.56, it could increase to 0.604 at night. The standard deviations of errors in 
daytime and nighttime surface temperatures, and band emissivities retrieved by the % 2 fit method are increased 
slightly. Note that the retrieved emissivities are compared to daytime emissivities only. In tests C4 and C5, we 
set different BRDF anisotropic factors for the three bands in the mid-infrared range by differences of 5% and 
10%. There is no significant change in the retrieved surface temperature and emissivities. In tests C6 and C7, we 
use non-Lambertian reflectance for the surface-reflected solar diffuse irradiance and atmospheric downward 
irradiance terms. They differ from the reflectance of a Lambertian surface by ±20%. The effect of the non- 
Lambertian reflectance is also not significant. Comparing the standard deviations in tests C2 through C7 to those 
in test Cl shows that the maximum difference in standard deviations of errors in retrieved surface temperatures is 
0.17 K and the maximum difference in standard deviations of errors in retrieved band emissivities is 0.005, they 
are comparable to or smaller than the effects caused by NEAT and calibration errors of the instrument. 
Therefore, we do not need to understand the three assumptions of surface optical properties described in section 
2.1 as strict constraints to the new LST algorithm. 
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4.4 Sensitivities to NEAT and Calibration Errors 

In the fourth simulation experiment, we keep the same atmospheric and surface temperature parameters as the 
first experiment, but change NEAT and calibration error values in a series of tests, as shown in Table V. The first 
column in the table indicates the test number. Seven NEAT values for seven bands used in the new LST 
algorithm are listed in the second column block, and a systematic calibration error for all bands in the third 
column. Standard deviations (8T 4 ) and maximum errors (A T s ) of the retrieved daytime and nighttime surface 
temperatures are given in columns 4-7. The standard deviations of errors in retrieved emissivities for MODIS 
bands 31 and 32 are given in the last two columns. Comparison between test D1 and test D2 indicates that the 
effect caused by a systematic calibration error of 0.5% is comparable to the effect of the given NEAT values. 
Test D3 indicates that doubling the NEAT values increases the standard deviation of retrieved daytime surface 
temperature by approximately 0.2 °K. Comparing tests D4 and D5 to test D2 indicates that errors in retrieved 
surface temperatures and band emissivities become larger as the calibration error increases. In order to achieve 
the 1 °K requirement for the LST accuracy and to retrieve band emissivities in MODIS bands 31 and 32 at an 
accuracy of the 0.01 level, the calibration error should be smaller than 1%. The new LST algorithm requires 
small NEAT (large signal-to-noise ratio) and a high consistent calibration accuracy for the seven bands used. The 
split-window SST and LST algorithms also need these requirements for MODIS bands 31 and 32. However, the 
new LST algorithm needs these requirements over a much wider spectral range. 

5. Several Issues Related to Day/night Algorithm 
5. 1 Day/night Registration 

It is natural to link the MODIS day/night LST algorithm to the apparent thermal inertia (ATI) algorithm of the 
Heat Capacity Mapping Mission (HCMM) sensor [NASA, 1980; Short and Stuart, 1982; Price, 1985; Vukovich, 
1984; Majumdar and Bhattacharya, 1990]. The HCMM program was the first NASA research effort directed 
mainly toward observations on the thermal state of Earth’s land surface from an unmanned satellite. The HCMM 
spacecraft operated between April 1978 and September 1980. MODIS is a much more advanced instrument than 
HCMM in terms of the thermal remote sensing capability except for the spatial resolution, as shown in Table VI. 
The nature of the analog telemetry system onboard the satellite results in some variability in the quality of the 
HCMM data and the data quality was also influenced by the quality of the receiving station recording [NASA, 
1980]. Because there was only one TIR channel on HCMM, it was not possible to make atmospheric and 
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emissivity corrections unless in-situ measurements of atmospheric profiles and surface emissivity were available. 
The ATI algorithm requires data in the 12-hour day and night coverage for estimate the day/night temperature 
difference. So cloud cover was a hindrance to the accomplishment of objectives [Short and Stuart, 1982]. The 
MODIS day/night LST algorithm does not require the 12-hour day and night coverage, it can use daytime and 
nighttime data collected in several days as long as the surface emissivity does not change significantly. 
Therefore, the chance to have a pair of daytime and nighttime data both in clear-sky conditions will be much 
large. It does not assume the same atmospheric conditions for daytime and nighttime because separate 
atmospheric variables are used for day and night. 

The uncertainty in day/night registration of MODIS data may be a major error source for the day/night LST 
algorithm. Numerical simulations have been made to evaluate the sensitivity of the day/night LST algorithm to 
the uncertainty in day/night registration. We assume that a vegetation component is mixed with another 
terrestrial material in the 80-sample database. It is assumed that the daytime proportion of the vegetation 
component in a mixed pixel is 0.5 and the nighttime proportion varies from 0.5 to 0.2 for simulating the mis- 
registration effect. The daytime canopy temperature of the vegetation component is given three values: the same 
surface temperature as for the another component, 4 °K warmer or cooler than the surface temperature of the 
another component. Its nighttime temperature is assumed to be equal to the surface temperature of the another 
component. The mid-latitude summer atmosphere is used in this simulation study. Note that the same NEAT 
values and the systematic calibration error 0.5% used in typical numerical experiments for the day/night LST 
algorithm are also used in this sensitivity study. The rms and maximum errors in surface temperatures and band 
emissivities retrieved by the y 2 fitting day/night LST algorithm are shown in Table VII. The band emissivities of 
the mixed pixel are calculated from the band emissivities in the database and the proportions. The average 
values of calculated daytime and nighttime band emissivities are used as the target values for retrieval. Results 
in Table VII indicate that the day/night LST algorithm still works well as long as the nighttime proportion of the 
vegetation component is not smaller than 0.35, differing from the daytime proportion by less than 30%. This 
corresponds to a mis-registration 15% or 30%, depending on how the vegetation component distributed in the 
mixed pixel. 

5.2 The effect of surface emissivity changes caused by dew 

Dew is a common global phenomenon. By model calculations and measurements of dew occurrence, Janssen 
and Romer [1991] showed that there were 220 dewy nights and 1600 dewy hours in the Netherlands in 1987. Six 
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years of dew observations in the Negev Desert, Israel [Zangvil, 1996] show that the total monthly dew hours is 
145 h per month from August to January, and 80 h per month from February to June and that the average duration 
of dew per dew night appears to follow closely the length of the night. Scherm and Bruggen [1993] used a 
previously verified and validated model to calculate the frequency and duration of dew occurrence with hourly 
weather data over a growing season (April-October) in different climate regions of California. Simulation results 
indicated that high dew occurrence (78-93% of the days) and long mean dew duration (8.7-9.3 per day) at coastal 
stations, and low occurrence (23% of the days) and short duration (0.9 h per day or 3.7 h per dew event) at the 
interior valley site. The dew duration is sensitive to changes in humidity and cloud cover at the coastal sites and 
to humidity and wind speed at the interior valley site. An empirical model to estimate the duration of dew 
periods by the use of dew point depression, wind speed, and relative humidity (RH) as inputs [Gleason et al., 
1994], has an accuracy better by only 5-10% than a simple RH > 90% model [Sutton et al., 1984], This indicates 
that the near-surface relative humidity is the primary factor in dew occurrence. Because the weather data are 
usually not available in remote areas, it is not easy to calculate dew occurrence with the dew models. The near- 
surface humidity in the MODIS atmospheric profile product may be used to reject the dew occurrence. 
Considering potential large errors in the temperature and water vapor in the surface boundary layer in the MODIS 
atmospheric profiles product, RH < 50% may be used as a conservative condition for nondew hours. Then we 
only need to detect night dew in our LST algorithm in relatively wet conditions (RH > 50%). When dew occurs 
at night the surface emissivity may be changed significantly if a land cover has a low emissivity in dry 
conditions, for example, in the desert environment. This significant change in surface emissivity should be 
shown in TIR data. If surface emissivity at night differs significantly from the value in daytime, the major link 
between day and night observations used in the day/night LST method will be lost. So we need to find some 
additional information to help the LST retrieval with only the multiband night-only observation. The empirical 
relationship between the minimum emissivity (e min ) and the spectral contrast which is described by the 
maximum-minimum difference (MMD) is such an information. This relationship is used in the ASTER TES 
algorithm [ATBD- AST-03, 1996]. Based on 80 laboratory reflectance spectra of rocks, soils, vegetation leaves, 
snow, ice and water [Salisbury and D’Aria, 1992] and their mixed values with three different portions of 
vegetation (totally 320 samples), two empirical relationships between minimum emissivity (e m j n ) and the spectral 
contrast described by MMD (maximum and minimum difference) have been calculated. One is for five MODIS 
bands (bands 20, 22, 23, 29 and 31), and another for the seven MODIS bands used in the day/night LST 
algorithm, as shown in Fig. 6. These two relationships have been added to our day/night algorithm as an option. 
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We also made a series of numerical simulations with the enhanced day/night LST algorithm (called MMD-MIN 
y 2 fitting LST method hereafter) in cases with and without night dew. Seven tests have been made to show the 
performance of three options of this enhanced LST algorithm: night-only, day-only, and day/night combination. 
In the current simulations we assume that the emissivity of a dew film has the same value of water surface 
emissivity. We will measure its spectral emissivity during the course of dew occurrence. 

Simulation results with the statistical day/night LST algorithm and the MMD-MIN y 2 fitting LST algorithm are 
shown in Fig. 7 and Table VIII. We used the similar parameters as used in the third numerical experiment and 
assume that each of 80 samples mixed with the vegetation sample (no. 51 in Table II) in a proportion of 25 
percent. Note that the spectral reflectance data of vegetation samples (no. 50 to 62) were obtained by 
measurements of single leaves. According to our recent measurements and simulations with modified BRDF 
kernel models [Snyder and Wan, 1996(a)], the spectral reflectance of a canopy is approximately one third of the 
spectral reflectance of a single leaf because of volumetric effects. We have made this adjustment for the 
vegetation samples in this numerical experiment in order to get a better understanding of the dew effect. After 
this adjustment, the emissivity of all vegetation samples is close to 0.99 so that the emissivity difference between 
water and vegetation is negligible. Figure 7 show the error in surface temperature and the emissivity in MODIS 
band 20 retrieved by the night-only method. For sample 50 to sample 79, the errors in nighttime surface 
temperature retrieved by both the statistical and y 2 fitting methods are within the range of ±1.5°K because dew 
does not change the surface emissivities of vegetation, water, and ice. But dew has a significant effect on the 
emissivities of most rock and soil samples as shown in Fig. 2A. As indicated in Table VIII, test A is for the 
surface temperature and emissivities retrieved by the day/night statistical method. Because the dew effect has 
not been included in the statistical regression analysis, the rms and maximum errors in nighttime surface 
temperatures are 1.7°K and 6.6 °K, respectively. In test B, only the night observation in seven MODIS bands 
was used to retrieve the surface temperature and emissivities by the y 2 fit method. The maximum error in the 
retrieved nighttime surface temperature reduces to 5.4 °K if no MMD-MIN relationship is used, and it reduces to 
3.4 °K and 3.1 °K if one or two MMD-MIN relationships are used. The rms error in retrieved band emissivities is 
reduced by a factor of 2-3. Figure 7A shows the error in surface temperatures retrieved by the statistical method 
in test A and by the y 2 fitting method in test B. Figure 7B shows the error in surface emissivities in band 20, £20 . 
retrieved by the the y 2 fitting method in test B. Although errors in the retrieved £20 are still large for some rock 
and soil samples, these errors are smaller than the emissivity difference between dew and these samples, which is 
shown by open circles in Fig. 7B. This means that £20 retrieved by the night-only method can be used to detect 
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dew for those land types which have low £20 values in dry conditions. Then we can use this information in the 
following processing. Test C shows the simulations results retrieved by the night-only x fitting LST method for 
cases without night dew, the accuracies are much better than those in Test B for cases with night dew. Tests D 
and E show the results retrieved by the day-only y} fitting method for cases with and without dew. The MMD- 
MIN relationships improve the retrieval accuracies slightly. The accuracies of retrieved emissivities in bands 20, 
22, and 23 are poor because the solar effects cannot be corrected well by the use of only the daytime observation 
data. In Test F for cases with night dew, use of both daytime and nighttime data does not improve the accuracies, 
even with the dew emissivity for nighttime. This is because the surface emissivity may be changed significantly 
between daytime and nighttime. Test G is the results retrieved by the day/night X fitting method with and 
without the MMD-MIN relationships in cases without night dew. Use of the MMD-MIN relationships in Tests F 
and G improves the accuracies by only a small amount. From the simulation results for Tests A through G, we 
can gain the following insights into the capability of the day/night LST algorithm to deal with the dew effects: 1) 
dew does not affect the LST and emissivity retrieval for most of the Earth surface which is covered by water and 
vegetation. 2) dew do affect the retrieval for land surfaces with low emissivities in dry conditions, mainly in 
semi-arid and arid areas. In these areas, dew reduces the accuracies in retrieved LST and emissivities in bands 
20, 22, and 23 roughly by a factor of 2, and reduces the emissivity accuracies in other bands slightly. 
Considering the low frequency of dew occurrence in semi-arid and arid regions, dew is not a serious problem for 
LST and emissivity retrieval but it makes the processing complicated slightly. We assume that the dew effect is 
needed to considered for one third of the global land (determined by the MODIS land-cover and NDVI products) 
at a frequency less than 50%. Because the dew detection process with the night-only LST mode takes less than 
half the computing time for the day /night LST mode, the total computing time will be increased by less than 10 
percent. 3) the MMD-MIN y} fitting day/night LST algorithm can be used to deal with night-only and day-only 
observations at lower accuracies. Surface emissivities in MODIS bands 20, 22, and 23 cannot be retrieved 
accurately by the day-only option in all situations (with and without dew). It is also worth noting the differences 
between the ASTER TES algorithm and the MODIS day/night LST algorithm: 1) the ASTER algorithm uses five 
bands in the spectral region 8- 12 pm, the MODIS algorithm uses four bands in the spectral range 8- 13.5 pm and 
three bands in the spectral region 3. 5-4. 2 pm. Two of the seven MODIS bands, bands 23 and 33, are on the edges 
of the atmospheric windows. They are used to adjust the boundary atmospheric temperature and the column 
water vapor but not for retrieval of surface temperature and emissivities per se. 2) If the atmospheric temperature 
and water vapor profiles recovered from the MODIS sounding channels are accurate enough (5T < 1 °K, the 
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accuracy of cwv < 10%), the ASTER algorithm can retrieve surface temperature and emissivity quite accurately. 
But this algorithm does not have the capability to reduce the uncertainties in the ASTER atmospheric correction 
product (including transmission, path radiance, and downwelling irradiance) which uses the MODIS atmospheric 
profiles product as input. The MODIS day/night LST algorithm uses the MODIS atmospheric profiles product 
only as initial conditions and has some capability to reduce the uncertainties in the atmospheric profiles by the 
use of self consistency in the seven MODIS bands. For water surfaces and evergreen canopies with stable and 
known emissivities, the LST and surface emissivity retrieved by this MODIS day/night algorithm will be 
compared with the values retrieved by the MODIS generalized split-window LST algorithm for additional 
evaluation of the self consistency. 3) The ASTER algorithm can retrieve daytime and nighttime surface 
temperatures and emissivities independently so that its performance will not be affected by the dew occurrence in 
principle. The dew occurrence will affect the MODIS day/night LST algorithm for some land-cover surfaces. 
Therefore, its night-only option will be used first for dew detection, and then compare the retrieved £20 to its 
daytime value retrieved previously. If there is a significant difference, this nighttime observation data should be 
avoided in the day/night algorithm. If this is the only nighttime observation available in one or two weeks, it can 
be used for day/night retrieval at a lower accuracy and it will be flagged in the QA (quality assessment) bits. 4) 
These two algorithm will produce LST and surface emissivities at different spatial resolutions. Comparison 
between the products by these two independent algorithms can provide useful information for validation. 

5.3 Cloud Cover 

Cloud cover is a common problem for visible and infrared remote sensing. Although the day/night LST 
algorithm can be used for all land surfaces, the main purpose of this algorithm is to retrieve the surface 
temperature and emissivity in semi-arid and arid regions where the surface emissivity varies in a wide range. The 
International Satellite Cloud Climatology Project (ISCCP) C2 data from 1984 to 1987 show that the world’s 
desert regions, which result from the sub-tropical anticyclones, are areas of total cloud amount minima, less than 
30% or 50% depending on the season [Drake, 1993]. From analysis of weather satellite observations, Chahine 
[1995] shows that the decreases in moisture and cloudiness are coupled with the increase in surface temperature 
over 304 °K, suggesting a positive feedback from the atmosphere perpetuating the existing desert conditions over 
dry and hot deserts. With EOS AM and PM platforms, the chance to have pairs of clear-sky day/night data will 
be increased significantly, especially in the dry areas. 
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6. Validation of the LST Algorithms 

We validated the physical principle of the day/night LST method with a ground-based sun-shadow method by 
measuring samples of soil, sands, grass and a black aluminum plate on the roof platform of our building at UCSB 
on January 19th and 26th, 1996. The solar beam is blocked for half of the samples. The TIR spectrometer views 
the portions in sunshine and in shadow for two separate measurements and also views a diffuse reflecting gold 
plate in the same spots for providing information of the solar and atmospheric downwelling radiation. After 
calibrating the spectrometer with blackbody at three different temperatures, another two separate measurements 
are made. For each sample, we obtained two pairs of data for the sunshine and shadow portions, and the diffuse 
reflecting gold plate. A band average procedure with the spectral response functions in seven MODIS TIR bands 
(i.e., 20, 22, 23, 29, 31-33) is used to achieve a high signal-to-noise ratio. Radiometric calibration is made with 
three blackbody temperatures, spectral emissivity of the blackbody surface and the front mirror, and the 
temperature of the front mirror. Then we use two methods to recover the surface temperature. In the 
conventional method, we use the spectral emissivity curves of samples measured with the integrating sphere 
system. In the sun-shadow method, we make non-linear % 2 fit of the sun-shadow data set for recovering surface 
temperatures in sunshine and in shadow, and the band-averaged emissivities. The LST values of samples of sand, 
soil, grass and black plate in sunshine and in shadow recovered by these two methods are shown in Fig. 8. Note 
that the mark squares represent the first method. The standard deviations are 0.4 °K and 0.1 °K, the maximum 
LST differences are 0.7 °K and 0.2 °K, for the LST difference in sunshine and in shadow, respectively. 

We conducted a field campaign jointly with the JPL (Jet Propulsion Laboratory) ASTER team at a large flat silt 
playa in Railroad Valley, Nevada, on August 3rd, 1995. MAS and TIMS (Thermal Imaging Multispectral 
Spectrometer) data, and field measurement data of surface spectral emissivity and temperature with TIR 
spectrometer and broadband radiometer were collected. Temporal and spatial analysis has been made. As shown 
in Table IX, LST retrieved from MAS data with the generalized split-window LST algorithm at view angle (0 V ) 
18.75° agrees with field measurement LST values within 1 °K. In this case, the LST accuracy is mainly limited 
by the uncertainty in its spatial variation. The MAS was calibrated with the new method [King et al., 1996]. 
Three field campaigns have been conducted in 1996. The first one was conducted over a snow field at the test 
site in Mammoth Lakes, California, on April 2nd, 1996. The second one was conducted jointly with other EOS 
teams at the same silt playa site in Railroad Valley, Nevada, on June 4th, 1996. The third one was conducted in 
the southern BOREAS test site on 14 August, 1996. MAS data and field measurement data were collected in 



28 


early afternoon and evening during the first two field campaigns. It was a perfect clear-sky day during the field 
campaign at the playa site on 4 June 1996. Surface temperatures were measured by three different methods: a 
TIR spectrometer, an infrared thermometer, and a thermistor 1mm beneath the surface. The surface temperature 
is also recovered from MAS data by the use of the generalized split-window and day/night LST algorithms. 
These results are shown in Fig. 9. Except the daytime temperature measured by the thermistor beneath the 
surface (it may differ from the surface temperature because of the variation in surface energy balance caused by 
the changing wind speed), the surface temperature values given by five methods are all within 1 °C. The band 
emissivities retrieved by the day/night LST method are lower than the values measured from playa samples, 
which were collected at the test sites, by 0.02 in bands 31-33, and 0.09 in bands 20, 22, and 23. Note the 
significant fluctuations in the daytime surface temperature caused by the change in wind speed. 


7. Computer Coding of the LST Algorithms 

Version 1 of the MODIS LST software has been submitted to the MODIS Team. 

In this version of the MODIS LST code, both the generalized split-window LST algorithm and the day/night LST 
algorithm have been implemented. The generalized split-window LST algorithm is used to retrieve LST for 
pixels with known band emissivities. The band emissivities are inferred from the land-cover type based on the 
previous month s MODIS land-cover product. The IGBP-type 1-km global vegetation database will be used at- 
launch as an alternative to the MODIS land-cover data. The numerical form of the generalized split-window LST 
algorithm is a linear equation of the band brightness temperatures of MODIS bands 31 and 32. A built-in look-up 
table in steps of 0.1 °K will be used to convert the band radiance into band brightness temperature. This LST 
algorithm will be very efficient in the processing. We also established a database of atmospheric look-up tables 
for MODIS bands 20, 22, 23, 29, 31-33, for the day/night LST algorithm. The data-base obtained from accurate 
atmospheric radiative transfer simulations and look-up tables and interpolation methods are used in these 
algorithms for numerical efficiency. The bubble chart for the LST generation process is shown in Fig. 10. The 
LST process is working on the granule basis. Each granule covers four or more tiles. The global sinusoidal grids 
system is divided into a large number of tiles, each with a nominal size of 1 100 km by 1 100 km (10° by 10° on 
the equator). An internal file of intermediate data, named MOD.AMl.Vl.lst_update.TILE, will be used in the 
process MOD_PRl 1 for each tile. This file stores the most updated clear-sky data for brightness temperatures in 
MODIS bands 20, 22, 23, 29, 31-33, atmospheric temperature and water vapor profiles at lower levels, date and 
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observation time, view angle, and related QA flags for both daytime and nighttime, and solar zenith angle, NDVI, 
condensed information of land-cover for daytime. The at-launch IGBP-type global vegetation database and 
NDVI will be used to distinguish pixels (class A) covered by dense vegetation, water, snow and ice from others 
(class B) in barren lands, including bare soils, semi-arid and arid regions. The generalized split-window LST 
algorithm will be applied to class A pixels for generating the level-2 LST product at 1km resolution. As doing an 
experiment, we used the global land-cover classification database (named VEG_CLSS) in the Global Data Sets 
for Land-Atmosphere Models, ISLSCP Initiative 1: 1987-1988, Volumes 1-5, which is available in CD-ROM 
from NASA Goddard DAAC Science Data Center. The spatial resolution of this database is 1 by 1 degree 
latitude/longitude. We can delineate semi-arid and arid regions which generally agree well with the dry regions 
defined by climatology [Trewartha and Horn, 1980]. 

The major steps in the LST process for each granule of 100 MODIS scan cube (1000 scan lines with 1354 1km- 
resolution samples each line) are described in the following. 

1) input one tile of the previous quarter’s land cover product and the intermediate file MOD.AMl.lst_update.TILE. 

2) input N scan-lines of MODIS geolocation I A and IB data and level-2 MODIS products (cloud mask , atmospheric 
profiles, snow , NDVI) used as input to the LST process. 

3) search for clear- sky pixels by the use of the MODIS cloud mask product, skip cloudy pixels. 

4) find the grids covered with each clear-sky pixel , and calculate their coverages. 

5) determine pixels with known band emissivities by the use of the previous quarter’s land-cover product, MODIS 
snow cover, and NDVI products. These products can be used to identify dense vegetation (evergreen needleleaf 
forest, evergreen broadleaf forest, deciduous needleleaf forest, deciduous broadleaf forest, grasslands, croplands, if 
NDVI > 0.8), water, snow, ice. Their emissivities in bands 31 and 32 can be calculated from the coverages of these 
land cover types and the emissivity knowledge base. Then LST can be retrieved with the generalized split-window 
LST algorithm. 

5’) for other pixels, retrieve LST with the generalized split-window algorithm by the use of the emissivity values in 
bands 31 and 32 that were retrieved by the day/night algorithm in previous days. This retrieved LST needs to be 
compared with and may be refined by the values retrieved by the day /night algorithm later. 

6) calculate and accumulate the values at each grid that will be used in the day/night LST algorithm. 

7) back to (2) if the granule is not completed. 

8) find new pairs of daytime and nighttime co-registered pixels 

9) proceed the statistical regression approach for the day/night LST method to retrieve surface temperatures and 
emissivities. 



30 


2 

10) proceed the X fi t approach of the day/night LST method for pixels in which the MODIS atmospheric 
temperature and water vapor profiles product has high quality. 

11) L3 LST quality assessment The LST retrieved by the day/night algorithm will be compared with and may be 
refined or replaced by the values retrieved by the generalized split-window algorithm, especially for land cover types 
with high and stable emissivities in bands 31 and 32. 

12) output L3 LST for the tile . 

13) update and output the intermediate tile file. 

14) L2 LST quality assessment. 

15) output N scan-lines ofL2 LST 

16) go to 1) if there is more tiles. 

The value of N in 2) is determined by the memory requirement and the size of the main memory on the CPU. We 
use 10 as its initial value. If N equals to the total number of scan lines in the granule, step 7) will be skipped. If 
the computer system has enough main memory, all tiles can be processed with only one-time input of the granule 
files. But most likely it is not the case in the beginning of the at-launch processing. 

8. Other Activities 

An International Land-Surface Temperature Workshop was held on September 17-19, 1996, at the University of 
California Santa Barbara. Jeff Dozier, Dean of the UCSB School of Environmental Science & Management 
welcomed the participants. Twenty five participants from USA, France, Australia and Japan attended the 
workshop. Twenty presentations were followed by two discussion sessions. It was a successful and productive 
workshop. The workshop was part of a continuing effort to maintain contact among members of the EOS 
community that are concerned with the improvement of LST algorithms, the definition of procedures for 
validation of LST, and the identification of the sources and the magnitude of measurement uncertainties. The 
specific goals of the Workshop were to clarify the present state of the art in LST estimation from spacebome 
sensors and to identify future directions including issues requiring further research effort. A subsidiary goal was 
to establish a closer relationship between LST algorithm designers and the LST user community. 

The new version (October 1996) of MODIS LST ATBD (algorithm theoretical basis document) was reviewed by 
five reviewers in November 1996 and by the EOS Review Panel on December 1 1, 1996. The scores given by the 
reviewers (2 As, 1 A-, and 2 Bs) reflect the significant improvement of the ATBD and the progress of the MODIS 
LST algorithms in the past two years. 
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Fig. 5, Histogram of errors in surface temperatures (A) and emissivities (B) retrieved by the X fit method. 
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TABLE I. Specifications of the EOS MODIS bands. 


band 

bandwidth 

(nm) 

IFOV 

primary 

use 

band 

bandwidth 

(pm) 

IFOV 

NEAT 

(°K) 

primary 

use 

1 

620- 670 

250m 

L 

20 

3.660-3.840 

1km 

0.05 

O, L 

2 

841-876 

250m 

A, L 

21 

3.929-3.989 

1km 


fire, volcano 

3 

459- 479 

500m 

L 

22 

3.929-3.989 

1km 

0.07 

A, L 

4 

545- 565 

500m 

L 

23 

4.020-4.080 

1km 

0.07 

A, L 

5 

1230-1250 

500m 

L 

1 

24 

4.433-4.498 

1km 

0.25 

A 

6 

1628-1652 

500m 

A, L 

25 

4.482-4.549 

1km 

0.25 

A 

7 

2105-2155 

500m 

A, L 

27 

6.535-6.895 

1km 

0.25 

A 

8 

405- 420 

1km 

O 

28 

7.175-7.475 

1km 

0.25 

A 

9 

438- 448 

1km 

O 

29 

8.400-8.700 

1km 

0.05 

L 

10 

483- 493 

1km 

O 

30 

9.580-9.880 

1km 

0.25 

ozone 

11 

526- 536 

1km 

O 

31 

10.780-11.280 

1km 

0.05 

A, L 

12 

546- 556 

1km 

O 

32 

11.770-12.270 

1km 

0.05 

A, L 

13 

662- 672 

1km 

O 

33 

13.185-13.485 

1km 

0.25 

A, L 

14 

673- 683 

1km 

O 

34 

13.485-13.785 

1km 

0.25 

A 

15 

743- 753 

1km 

A 

35 

13.785-14.085 

1km 

0.25 

A 

16 

862- 877 

1km 

A 

36 

14.085-14.385 

1km 

0.35 

A 

17 

890- 920 

1km 

A 






18 

931-941 

1km 

A 






19 

915- 965 

1km 

A 






26 

1360-1390 

1km 

cirrus 







Note: A - atmospheric studies; L - land studies; O - ocean studies. Ref: MODIS Level IB Algorithm Theoretical 
Basis Document, 1995, NASA/GSFC, Greenbelt, MD. 
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TABLE II. List of terrestrial material samples. 


sample 

no. 

sample 

name 

type of 
material 

sample 

no. 

sample 

name 

type of 
material 

1 

basalt, f 

fresh rough surface 

41 

0135 

soil (Entisols) 

2 

basalt.v 

desert vanish coated rock 

42 

0145 

soil (Ultisols) 

3 

ijolite.f 

fresh rough surface 

43 

0211 

soil (Molisols) 

4 

ijolite.v 

desert vanish coated rock 

44 

0219 

soil (Alfisols) 

5 

rhyolite.f 

fresh rough surface 

45 

0226 

soil (Inceptisols) 

6 

rhyolite.v 

desert vanish coated rock 

46 

0475 

soil (Vertisols) 

7 

crustose.10 

lichens coated rock 

47 

1530 

soil (Aridisols) 

8 

crustose.65 

lichens coated rock 

48 

4717 

soil (Oxisols) 

9 

basalt.h7 

igneous rock 

49 

foliose.l 

veg., lichens 

10 

dunite.hl 

igneous rock 

50 

indiangr.ass 

veg., green foliage 

11 

granite. hi 

igneous rock 

51 

redoak 

veg., green foliage 

12 

syenite.hl 

igneous rock 

52 

white, ine 

veg., green foliage 

13 

greywack.ehl 

sedimentary rock 

53 

senbeech 

veg., senescent foliage 

14 

limeston.ehl 

sedimentary rock 

54 

senpine 

veg., senescent foliage 

15 

limeston.eh2 

sedimentary rock 

55 

senredoa.khl 

veg., senescent foliage 

16 

limeston.eh3 

sedimentary rock 

56 

senryegr.ass 

veg., senescent foliage 

17 

sandton.ehl 

sedimentary rock 

57 

oakbark.l 

veg., tree bark 

18 

sandton.eh2 

sedimentary rock 

58 

pinebark.1 

veg., tree bark 

19 

sandton.eh4 

sedimentary rock 

59 

ypoplarb.ark 

veg., senescent foliage 

20 

shale.h3 

sedimentary rock 

60 

coniferous 

veg. decomposing litter 

21 

shale.h5 

sedimentary rock 

61 

decidu.ous 

veg. decomposing litter 

22 

shale.h6 

sedimentary rock 

62 

wood 

veg. decomposing litter 

23 

siltston.ehl 

sedimentary rock 

63 

seawater 

water 

24 

siltston.eh2 

sedimentary rock 

64 

distwa.ter 

water 

25 

gneiss. hla 

metamorphic rock 

65 

disticel .00g 

ice 

26 

gneiss. h3a 

metamorphic rock 

66 

distices.moo 

ice 

27 

gneiss. h4 

metamorphic rock 

67 

seaice.10.ogr 

ice 

28 

marble. h2 

metamorphic rock 

68 

seaicesm.oot 

ice 

29 

marble. h3 

metamorphic rock 

69 

qtzwater.23 

suspended sediments 

30 

marble. h4 

metamorphic rock 

70 

qtzwater.64 

suspended sediments 

31 

quartzit.ehl 

metamorphic rock 

71 

qtzwater.7 

suspended sediments 

32 

quartzit.eh4 

metamorphic rock 

72 

foam 

water coatings 

33 

quartzit.eh6 

metamorphic rock 

73 

oil15465 

water coatings 

34 

schist.h3a 

metamorphic rock 

74 

oil34792 

water coatings 

35 

schist.h6a 

metamorphic rock 

75 

oil39076 

water coatings 

36 

schist.h7 

metamorphic rock 

76 

oil42667 

water coatings 

37 

slate. hi a 

metamorphic rock 

77 

soilfl.oat 

water coatings 

38 

slate. h2a 

metamorphic rock 

78 

qtzfloat 

water coatings 

39 

slate.h3 

metamorphic rock 

79 

oil35473 

water coatings 

40 

0127 

soil (Spodosols) 

80 

qtz.hem 

quartz 
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TABLE III. The standard deviations of errors in surface temperature and emissivities retrieved with two 
approaches of the day/night LST algorithm. 


atmospheric 

profile 

8 Tsday 
(° K) 

8 Ts- night 

(°K) 

8 £20 

hZ22 

8e23 

8E29 

&31 

8 c 32 

8C33 



with the statistical regression approach 




average 

0.91 

0.73 

0.021 

0.025 

0.027 

0.013 

0.012 

0.014 

0.012 

A109 

0.82 

0.75 

0.026 

0.024 

0.027 

0.032 

0.013 

0.018 

0.014 

A115 

1.18 

0.73 

0.021 

0.027 

0.033 

0.015 

0.013 

0.016 

0.013 

average-4K 

0.94 

0.64 

0.019 

0.022 

0.024 

0.013 

0.013 

0.014 

0.013 




p 

with the X fit approach 





average 

0.51 

0.36 

0.015 

0.014 

0.016 

0.008 

0.009 

0.009 

0.012 

A109 

2.13 

1.91 

0.057 

0.068 

0.076 

0.030 

0.036 

0.043 

0.014 

A115 

0.97 

0.58 

0.028 

0.024 

0.032 

0.023 

0.017 

0.020 

0.013 

average-4K 

0.45 

0.32 

0.014 

0.013 

0.015 

0.008 

0.009 

0.009 

0.013 


TABLE IV. The RMS errors in surface temperature and emissivities retrieved with the x 2 fit approach of 
the day/night LST algorithm in the sensitivity study on assumptions of surface optical 
properties in conditions of T a - ( j a y=298.2°K, T a-night =290.2°K, cwv = 2.6cm, a = 1.0, 
NEAT =0.05-0. 12° K, and systematic calibration error = 0.5%. 


test 

test 

8 T s day 

8 Ts- night 

§£31 

8 £32 

no. 

conditions 

(°K) 

(°K) 

(daytime) 


£n(j)=£d(j) 





Cl 

1 

li 

II 

0.51 

0.36 

0.009 

0.009 


Lambertian surface 





C2 

Zn(i)=Zd(j) + 0.01 

0.75 

0.55 

0.013 

0.013 

C3 

£n(j)=£d(j) + 0.1(1 -Ed) 

0.51 

0.41 

0.009 

0.009 

C4 

f r (22) = 0. 95 f r (23), f r (20) = 0. 90 f r (23) 

0.71 

0.58 

0.012 

0.011 

C5 

f r (20)= 1. 10f r (23), f r (22)= 1.05f r (23) 

0.68 

0.60 

0.011 

0.011 

C6 

non-Lambertian surface ( 80%) 

0.61 

0.43 

0.011 

0.010 

C7 

non-Lambertian surface (120%) 

0.58 

0.41 

0.010 

0.010 
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TABLE V. The dependences of standard deviations (5T S ) and maximum errors (A T s ) in surface 
emissivities and temperatures retrieved with the % 2 fit approach of the day/night LST 
algorithm on NEAT and calibration errors. 


calibration § T s _ d a y 8 T~s- night A T s (j a y A 7~ s _ n; g^f $£ 3 / &32 


no. 

( °K) 

errors (%) 

(°K) 

(°K) 

(°K) 

(°K) 



D1 

0.05.0.07,0.07.0.05,0.05,0.05,0. 1 2 

0.00 

0.41 

0.31 

3.3 

2.6 

0.007 

0.007 

D2 

0.05,0.07,0.07,0.05,0.05,0.05,0. 1 2 

0.50 

0.51 

0.36 

3.2 

2.1 

0.009 

0.009 

D3 

0.10,0.14,0.14,0.10,0.10,0.10,0.25 

0.50 

0.69 

0.49 

3.7 

2.2 

0.011 

0.012 

D4 

0.05,0.07,0.07,0.05,0.05,0.05,0. 1 2 

0.75 

0.58 

0.40 

3.3 

2.1 

0.010 

0.011 

D5 

0.05,0.07,0.07,0.05,0.05,0.05,0.12 

1.00 

0.66 

0.45 

4.4 

2.5 

0.012 

0.012 


TABLE VI. Comparison between HCMM and MODIS TIR channels. 


specification 

HCMM 

MODIS 

orbital altitude 

620km 

705km 

resolution 

0.6km at nadir 

1km at nadir 

swath width 

716km 

2330km 

spectral range 

10.5-12.5pm 

3.5-14.5|im 

number of TIR channels 

1 

16 

NEAT 

0.4 °K 

0.05 °K for surface channels 



0.25 °K for sounding channels 

quantization 

8 bits from analog data 

12 bits 

positioning error 

= 3km 

< 200m 

spatial coverage 

direct broadcast to ground stations 

gloabl 

launch date 

April 1978 

June 1998 

life time 

17 months 

5 years 
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TABLE VII. The sensitivity of surface emissivities and temperature retrieved by the x 2 fit day/night LST algorithm 
to uncertainties in day/night registration. The proportion of vegetation within a mixed pixel in daytime 
differs from that in nighttime because of mis-registration. 


P(veg) 


rms errors 


maximum errors 


day 

night 

87 * s-ct 

(°K) 

8Tg_ n 

( °K) 

5£20 

&3i 

5e 32 

ATg-d 

( °K) 

ATg_ n 

(°K) 

At20 

Ae 31 

Az 32 



T s~d( veg) - Ts_tf (' background ) =-4°K 







0.5 

0.5 

0.35 

0.24 

0.011 

0.007 

0.008 

3.1 

1.5 

0.090 

0.041 

0.040 

0.5 

0.45 

0.35 

0.25 

0.012 

0.008 

0.009 

1.8 

1.2 

0.059 

0.029 

0.035 

0.5 

0.4 

0.37 

0.29 

0.014 

0.008 

0.010 

1.8 

2.2 

0.092 

0.035 

0.041 

0.5 

0.35 

0.44 

0.36 

0.016 

0.010 

0.012 

2.6 

3.4 

0.130 

0.060 

0.060 

0.5 

0.3 

0.46 

0.41 

0.018 

0.011 

0.014 

2.9 

3.5 

0.145 

0.074 

0.078 

0.5 

0.25 

0.56 

0.51 

0.022 

0.013 

0.017 

4.2 

4.5 

0.182 

0.086 

0.081 

0.5 

0.2 

0.61 

0.57 

0.023 

0.015 

0.018 

4.1 

3.4 

0.185 

0.144 

0.101 



T^f veg) - Tg_ ^ (' background ) = 0°K 







0.5 

0.5 

0.36 

0.26 

0.012 

0.007 

0.009 

2.3 

1.5 

0.068 

0.031 

0.037 

0.5 

0.45 

0.38 

0.27 

0.013 

0.008 

0.010 

2.1 

1.5 

0.068 

0.037 

0.037 

0.5 

0.4 

0.38 

0.30 

0.014 

0.009 

0.011 

1.8 

1.9 

0.082 

0.036 

0.038 

0.5 

0.35 

0.43 

0.36 

0.016 

0.010 

0.013 

2.6 

2.8 

0.109 

0.047 

0.048 

0.5 

0.3 

0.46 

0.44 

0.019 

0.012 

0.015 

2.9 

4.7 

0.183 

0.083 

0.074 

0.5 

0.25 

0.52 

0.51 

0.021 

0.013 

0.017 

3.7 

3.6 

0.171 

0.068 

0.077 

0.5 

0.2 

0.58 

0.59 

0.023 

0.015 

0.019 

4.4 

4.3 

0.150 

0.185 

0.091 



Tg-df veg) - Tg^ (' background ) = 4° K 







0.5 

0.5 

0.39 

0.25 

0.011 

0.007 

0.009 

3.2 

1.5 

0.088 

0.037 

0.037 

0.5 

0.45 

0.40 

0.26 

0.012 

0.008 

0.010 

2.4 

1.4 

0.069 

0.033 

0.037 

0.5 

0.4 

0.45 

0.32 

0.014 

0.009 

0.012 

2.6 

2.3 

0.096 

0.039 

0.046 

0.5 

0.35 

0.50 

0.38 

0.016 

0.011 

0.014 

2.9 

3.3 

0.128 

0.059 

0.059 

0.5 

0.3 

0.57 

0.45 

0.019 

0.012 

0.016 

4.0 

4.5 

0.178 

0.083 

0.068 

0.5 

0.25 

0.65 

0.54 

0.021 

0.014 

0.018 

4.7 

4.3 

0.163 

0.072 

0.062 

0.5 

0.2 

0.73 

0.62 

0.024 

0.015 

0.020 

5.3 

3.9 

0.182 

0.063 

0.068 
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TABLE VIII. The standard deviation (5) and maximum (A) errors in surface temperature and emissivities retrieved 
with the statistical and the MMD-MIN x 2 -fitting LST algorithms in cases with and without night dew. 


number of 
MMD-MIN 
relations used 

5 7 ~$-<jay 
(A 

1 

5 Ts-night &20 

A T$-night A &20 

(°K) 

22 

Ae.22 

8 E 23 

A£23 

5£29 

Ae 29 

&31 

&31 

bt 32 

Ae 32 

8 £33 
^33) 


day/night statistical LST 

method for cases with night dew (including daytime e’s only) 


0 

0.87 

1.71 

0.026 

0.029 

0.026 

0.023 

0.012 

0.013 

0.014 


(3.63 

6.56 

0.140 

0.132 

0.097 

0.113 

0.045 

0.052 

0.040) 


night-only x 

^-fitting LST 

method for cases with unknown night dew (retrieving dew e’s) 


0 


1.49 

0.070 

0.066 

0.071 

0.037 

0.033 

0.033 

0.007 



(5.37 

0.221 

0.204 

0.214 

0.123 

0.109 

0.116 

0.037) 

1 


0.84 

0.039 

0.038 

0.042 

0.019 

0.015 

0.014 

0.007 



(3.36 

0.161 

0.147 

0.149 

0.094 

0.062 

0.087 

0.037) 

2 


0.81 

0.037 

0.037 

0.043 

0.017 

0.013 

0.011 

0.007 



(3.14 

0.149 

0.137 

0.146 

0.094 

0.068 

0.071 

0.037) 


night-only x 

-fitting LST 

method for cases without night dew (retrieving nighttime e’s) 


0 


0.46 

0.020 

0.019 

0.021 

0.011 

0.009 

0.009 

0.012 



(1.75 

0.063 

0.058 

0.064 

0.057 

0.037 

0.034 

0.031) 

1 


0.44 

0.020 

0.019 

0.021 

0.012 

0.010 

0.011 

0.012 



(1.77 

0.058 

0.059 

0.063 

0.084 

0.039 

0.035 

0.031) 

2 


0.42 

0.020 

0.019 

0.021 

0.011 

0.009 

0.008 

0.012 



(1.81 

0.068 

0.075 

0.080 

0.092 

0.035 

0.037 

0.031) 


day-only x 2 

-fitting LST method for cases with unknown night dew (retrieving daytime e’s) 


0 

0.74 


0.021 

0.055 

0.059 

0.014 

0.012 

0.013 

0.014 


(2.74 


0.110 

0.304 

0.306 

0.055 

0.047 

0.045 

0.040) 

1 

0.63 


0.021 

0.028 

0.048 

0.011 

0.009 

0.010 

0.014 


(3.19 


0.111 

0.215 

0.296 

0.046 

0.041 

0.053 

0.040) 

2 

0.64 


0.021 

0.028 

0.047 

0.011 

0.009 

0.009 

0.014 


(2.69 


0.117 

0.210 

0.277 

0.041 

0.040 

0.043 

0.040) 


day-only % 2 - 

fitting LST method for cases without 

night dew 

(retrieving 

daytime e’s) 



0 

0.74 


0.015 

0.066 

0.057 

0.014 

0.012 

0.012 

0.012 


(2.42 


0.059 

0.314 

0.191 

0.074 

0.067 

0.063 

0.031) 

1 

0.57 


0.017 

0.035 

0.039 

0.009 

0.007 

0.008 

0.012 


(2.42 


0.291 

0.300 

0.265 

0.032 

0.037 

0.064 

0.031) 

2 

0.56 


0.016 

0.029 

0.040 

0.009 

0.007 

0.007 

0.012 


(2.42 


0.148 

0.287 

0.295 

0.044 

0.032 

0.052 

0.031) 


day/night % 2 

-fitting LST method for cases with known night dew (retrieving daytime e’s) 


0 

0.74 

0.21 

0.021 

0.055 

0.059 

0.014 

0.012 

0.013 

0.014 


(2.74 

2.09 

0.110 

0.304 

0.306 

0.055 

0.047 

0.045 

0.040) 

1 

0.63 

0.20 

0.021 

0.028 

0.049 

0.011 

0.009 

0.010 

0.014 


(3.19 

2.09 

0.111 

0.146 

0.296 

0.046 

0.041 

0.053 

0.040) 

2 

0.63 

0.20 

0.021 

0.030 

0.049 

0.011 

0.009 

0.009 

0.014 


(2.69 

2.09 

0.117 

0.255 

0.277 

0.041 

0.040 

0.043 

0.040) 


day/night x 2 - 

-fitting LST method for cases without 

night dew 

(retrieving 

averaged e’s) 



0 

0.40 

0.30 

0.014 

0.013 

0.014 

0.008 

0.007 

0.007 

0.012 


( 2.02 

1.49 

0.073 

0.065 

0.072 

0.110 

0.033 

0.034 

0.031) 

1 

0.37 

0.28 

0.013 

0.012 

0.013 

0.007 

0.006 

0.007 

0.012 


(1.75 

1.44 

0.054 

0.051 

0.057 

0.075 

0.029 

0.033 

0.031) 

2 

0.33 

0.26 

0.011 

0.011 

0.012 

0.006 

0.005 

0.006 

0.012 


(1.77 

1.20 

0.051 

0.055 

0.060 

0.053 

0.027 

0.031 

0.031) 
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TABLE IX. Summary of LST values over the test site (38’ 31.46’N, 115’ 42.74’W) in Railroad Valley, Nevada, 
during 1:22 and 1:30 PDT on 8/3/95. The size of one MAS pixel is approximately 50m by 50m. 


size of area 

mean (°C) 

stdv (°C) 

remarks 

12 cm diameter 

58.5 


by radiometer 

5 cm diameter 

59.2 


by spectrometer at 0^ 20 ° 

1 MAS pixel 

59.1 


ate„ 18.75° 

3 by 3 MAS pixels 

58.9 

0.48 


5 by 5 MAS pixels 

58.8 

0.67 


7 by 7 MAS pixels 

58.9 

0.76 


9 by 9 MAS pixels 

59.0 

0.81 


11 by 1 1 MAS pixels 

58.9 

0.82 


21 by 21 MAS pixels 

58.9 

1.21 




